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A B()i;rZMANN  TRANSPORT  CODE 
FOR  ION  PENETRATION  IN  MATTER 

A.  Introduction 

Boltzmann  transport  theory  has  been  thoroughly  investigated  for 
cases  where  specific  physical  assumptions  - isotropic  elastic  scattering 
cross  sections  in  the  center  of  mass  frame  (slow  neutron  transport)  or 
small  angle  elastic  scattering  (passage  of  high  energy  particles 
through  matter)  - allow  drastic  mathematical  simplifications.  No  such 
simplifying  conditions  exist  when  one  is  considering  either  low  energy 
ions,  which  come  to  rest  in  the  bombarded  material,  or  the  recoils  they 
produce  in  elastic  collisions  with  target  atoms.  However,  the  single 
particle  distribution  functions  for  the  penetrating  ions  and  the  recoil 
atoms  contain  a great  deal  of  information,  which  may  be  used,  e.g.,  to 
calculate  directly  range  distributions  of  both  beam  and  recoil  atoms, 
energy  and  damage  distributions  in  the  bombarded  material  and  to  study 
ion  implantation  and  sputtering  phenomena.  1.  Manning  and  D.  W.  Padgett^ 
(hereafter  referred  to  as  MP)  have  developed  a formalism  for  describing 
the  penetration  of  amorphous  matter  by  a heavy  ion  beam  which  is  based 
on  the  Boltzmann  transport  equation  and  uses  Lindhard  atomic  collision 
cross  sections.  This  formalism  was  later  extended  by  Mueller^  to 

Include  inelastic  losses. 

Notr:  Manuscript  submitted  July  6,  1977. 
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In  this  paper,  a method  in  wliich  the  single  particle  distribution 
function  is  expanded  in  Hermite  and  Legendre  polynomials  is  described 
and  applied  to  the  implementation  of  the  Manning-Padgett  heavy  ion 
transport  formalism.  By  way  of  a demonstration  calculation  we  have 
considered  a beam  of  2Cf  keV  antimony  ions  incident  on  germanium.  This 
particular  case  was  chosen  because  it  seems  typical  for  a class  of 
applications  we  have  in  mind,  and  because  experimental  data’’  and  theo- 
retical calculations^  of  energy  deposition  are  already  available  here. 

As  discussed  in  MP,  most  of  the  calculations  of  ion  penetration 
have  heretofore  been  based  on  transport  theories  of  the  Lindhard 
type,^’*^’  whereas  the  transport  theory  we  deal  with  is  based  on  the 
Boltzmann  transport  equation. 


B.  Form  of  the  Solution 

Consider  the  Boltzmann  equation  for  the  vector  flux  t of  the  inci- 
dent beam,  under  the  simplifying  assumptions  of  rotational  invariance, 
translational  invariance,  plane  symmetry  and  time  independence:*^ 


^ ■^(E,x,|u)  + Na'^(E)'i(E,x,u)  - (S(E)Y  (E,x,a)) 

= + Na  V J d-v'(  y'v')  a(v'  • vl  Y(E x,  , 


(1> 


where  x is  the  distance  along  the  direction  of  the  bombarding  beam  and 
perpendicular  to  the  target  surface,  v is  the  ion  velocity,  N is  the 
density  of  target  atoms. 
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E = g -nv'  (2) 

U = V • X . (3) 


S(E)  is  the  inelastic  stopping  power  (we  take  inelastic  electronic 
energy  losses  into  account  in  the  continuous  slowing  down  approximation), 
> is  tlie  source  term  of  the  beam  ion  flux,  and 


'f(E,x,^)  = uvN(E,x,u)  , (■'.) 

with 


N(  E,  X,  ^)  d E d X d = number  of  beam  ions  with 

energy  in  the  interval 

dE  about  E,  direction 

cosine  du  about  and 

at  distance  dx  about  x;  (5) 


and  the  elastic  binary  collision  cross  section  is  defined  such  that 


> t)  d"^v^  d^r  dt 


= the  probable  number  of  beam  particles  with 
velocity  in  d^v ' about  v%  located  in  d^r 
about  which  undergo  a collision  in  time 
dt  about  t such  that  their  final  velocities 
lie  in  d'v  about  (6) 

T 

The  total  elastic  cross  section,  a (E),  is  given  by 

o^(E)  = I d^v'o( V - v')  . (Y) 

It  should  be  stressed  that  our  transport  approach  models  the  target 
as  being  amorphous. 
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We  employ  the  boundary  condition  used  by  Winterbon,  Sigmund  and 
Sanders  (WSSl  in  tlieir  work^  and  model  the  target  as  being  infinite  in 
every  direction;  over  the  entire  x = 0 plane  there  is  embedded  in  the 
target  a source  of  monoenerget ic  atoms,  all  traveling  in  the  positive  x 
direction.  This  boundary  condition  represents  an  actual  ion  bombardment 
only  to  the  extent  that  scatterings  back  and  forth  through  the  x = C 
plane  can  be  neglected.  For  the  case  of  our  demonstration  calculation, 
there  should  be  a negligible  number  of  such  scatterings. 

We  note  however  that  our  approach  nevertheless  differs  to  a con- 
siderable extent  from  that  of  WSS.  The  WSS  distribution  F(r,vl  is  a 
distribution  of  deposited  energy,  whereas  ours  is  the  distribution  in 
energy  and  angle  of  the  beam  ions  as  they  traverse  the  target.  In 
addition,  WSS  solve  for  moments  of  their  distribution  from  which  they 
must  reconstruct  the  distribution,  whereas  we  solve  directly  for  the 
distribution  by  an  expansion  in  Hermite  and  Legendre  polynomials  with 
energy  dependent  parameters  built  into  the  expansion  functions.  The 
use  of  these  energy  dependent  functions,  x(E)  and  5(E)  (see  next 
section),  allows  us  great  flexibility  in  improving  the  convergence  of 
the  expansion.  Given  Y(E,x,|j),  we  are  able  to  calculate  with  relative 
ease  the  distribution  of  deposited  energy,  the  energy  dependent  spatial 
distribution  of  beam  ions  and  moments  of  these  distributions. 

A method  of  handling  the  singularities  near  the  boundary  surface 
is  presented  in  references  1 and  2.  As  is  done  there  the  solution  for 
the  vector  flux  is  written  in  the  form 

= w + Y Y < $ . 
o 1 p 

)* 


Y 


(3) 


The  terms  and  represent  the  vector  flux  of  particles  which 

have  undergone,  respectively,  zero,  one  and  two  elastic  collisions  with 
target  atoms.  The  ^ term  represents  particles  which  have  undergone 
three  or  more  elastic  collisions. 

All  of  our  results  are  normalized  to  unit  incoming  flux.  We  write 
the  source  term  in  the  Boltzmann  equation  for  the  incoming  ions  in  the 
form 


y(E,x,u)  = 46(x)6(u  - 1)6(E  - Eg)  , (9) 

.e  E_  is  the  incident  beam  energy. 

D 

For  convenience  we  rewrite  Eq.  (8)  so  as  to  exhibit  explicitly  the 
singularities  arising  from  the  source  term. 

Y(E,x,u)  = 6(u  - 1)6(E  - E„),^(x) 

D O 

+ 6(m  - g(Eg,E))y^(E,x) 

+ Y2(E,X,u) 

+ $(E,x,^)  , (10) 

where  g(E  ,E)  is  the  cosine  of  the  scattering  angle  in  the  laboratory 
f5 

system  for  a beam  ion  entering  with  energy  E and  exiting  with  energy  E. 

D 

The  first  three  contributions  are  solved  for  explicitly,  and  the  problem 
is  reduced  to  solving  for  i,  which  is  by  far  tlie  most  important  portion 
of  Y>  hut  contains  neither  the  singularities  of  the  Yq  and  Yj^  terms  nor 
the  extremely  peaked  behavior  that  Y^  exhibits.'^ 
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C.  Explicit  Expression',  lor 


. nnd  Y 


It  is  convenient  to  define  a number  of  quantities  here  which  appear 
frequently  in  many  expressions  that  follow.  For  an  elastic  scattering, 
corresponding  to  v'  • ;v  (wliore  both  v'  and  v refer  to  the  beam  ion), 
one  obtains  for  the  cosine  of  the  scattering  angle  in  the  lab  system 


O'-  0 = g(E',E)  = :[(M  t i)(E/E')"  - (M  - 1)(e'/E)'] 


where 


M = A /A 

d 1 


and  A-  and  Ag  ate  the  masses  of  the  beam  ions  and  target  atoms, 
respectively.  The  cross  section  a(v ' - of  Eq.  (6)  can  be  written 
in  the  form 


cr(v'  V)  = -^6(0-0'  - g(E',E))F(E',E) 


For  convenience,  we  also  define  the  quantity 


3I(e',E)  = (2/A^3)^  NE^  F(E',E)  . 


The  total  cross  section  corresponding  to  the  differential  cross  section 


(13)  is 


T 

No  (E')  = f dEtf(E',E)  , 
E'/(5 


where 


/ M -t  1 y 

= V M - 1 ; 
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ami  e'/  is  the  minimum  possible  energy  for  an  outgoing  beam  ion  follow- 
ing an  elastic  collision.  For  later  use,  we  also  define  the  quantities 

ECE^a)  = N -^(El/u  . (17) 

t(E)  = L(E,g(E^^,E))  = No'^(E)/g(E^,E)  , (iS) 

f ,5  = = Na'^’(Ej^)  . (13) 

In  evaluating  and  V.  we  have  neglected  the  contribution 

from  inelastic  electronic  energy  losses.  As  will  be  seen  later  (Fig.  I") 
this  should  have  a negligible  effect  on  the  final  solution.  Inelastic 
losses  were  not  neglected  in  the  calculation  of  $.  In  deriving  the 
equation  below  we  assume  the  inelastic  loss  term  in  the  Boltzmann 
equation  to  be  absent.  For  the  sake  of  clarity  this  is  true  also  for 
the  equations  that  determine  i.  In  Appendix  A,  however,  we  show  how 
tliose  equations  are  modified  to  include  the  contribution  of  these  losses. 

Wlien  we  substitute  the  solution  (1C)  into  the  Boltzmann  equation 
and  use  the  limiting  procedures  of  to  separate  the  various  terms 

according  to  their  singular  behavior,  we  obtain  the  four  equations 

,^(x)  + No^(E,^),^(x)  = 6(x)  , (20) 

g(E|j,E) V ,(E,x)  ♦ No’^(E)t^(E,x)  (21) 

= g(Eg,E)t^(x)'/(E,^,E)  , 

- ^ ■*  Nc^(E)Y^(E,x,h)  (22) 

,,„PE  ... 

- -T—  f dE  '‘.F(e',E)  f dv  '6(v  . v'  - g(E',El)  X 
c:tt  „ V 
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X u'"^6( 


/ 


- g(E„,E'))y  (E',x^; 

O 1 


a1  so 


^ ^(E,x,^)  t No’^(E)$(E,x,n)  ( .2) 

o ^ 

= dE'JCE'jE)  r dv'6(v  • v'  - g(E',E))  X 

E 

X la'■‘[Y2(E',x,^')  + $(E',x,^')J  . 


Equation  (20)  is  simply  solved  to  yield 


r,  (x)  = U(x)  exp  {-  x]  , 


(2'’0 


where  L'(x)  is  the  unit  step  function.  In  solving  for  ^ , we  note  that 
the  mass  ratio  N in  the  case  of  antimony  on  germanium  is  such  that  the 
maximum  scattering  angle  is  approximately  o6  . The  cosine  factor 
g(E  ,E)  in  Eq.  (21)  can  never  be  negative;  in  fact^  three  or  more  col- 

D 

lisions  are  required  before  the  kinematics  allow  a beam  ion  to  be  scat- 
tered to  a backward  direction. 

As  can  be  easily  verified  by  substitution,  the  solution  of  Eq.  (21) 
has  the  fonn 


t (E,x)  = e dx'c”^^^^^  y (x')?(E  ,E)  , 

1 c l> 


(25) 


or 


t^(E,x)  = 3:(Eg,E)T,(l.p,f,(E)  ;x)  , 


(26) 


where 
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-ax  -bx 

T (a,b;x)  = — — U(x)  . 

D - a 


(^7) 


Mass  ratios  that  allow  backward  scattering  in  a single  collision  are 
handled  just  as  simply,  but  the  resulting  expression  for  ^ is  slightly 
d if  feren t . 

By  performing  the  angular  integrals,  Eq.  (22)  can  be  rewritten  as 


Y^(E,x,  ^)  + Na’'’(E)Y  (E,x,  u) 


dE'5(E',E)^^(E',x)g-HEg,E')K‘2U(K)  , 
E 


.1 


(28) 


where 


< = k(e',e,u)  = (h  - u)(h  - U ) 


(29) 


and 


= g(Eg,E')g(E',E)  ± /rr^'(Eg,E')  /r-'g2(E',E) 


(3C) 


As  can  again  be  verified  by  substitution,  the  solution  of  Eq.  (29)  is 


V (E,x,n)  = ^ d E '3=(  E ' , E) 


n 


.1 


(31) 


X t^(E',x')g-l(Eg,E')K"^  UOO 


= 4 j‘^^dE'y(E%E)5(Eg,E')g-l(Eg,E') 
E 


(32) 


X K ,i(E'),  L(E,|i)  ;x)  , 


where 
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(I- 


T.  (a,b,c;x)  = 


I -ax  -cx 
' e - e 


-bx  -cx 
e - e 


b - a / c - a 


c - b \ ‘ 


After  some  manipulation,  one  can  show  that  K.(e',E,,^)  can  be  written 
in  the  form 


(E  - e'Ue'  - e ) ‘ (M‘  - l)n 

K(E',E,a)  = ' ° 


e' 


?/  EE 


(3i*) 


B 


where 


U = + 1)^(E/E^)^  + (M  - 1)2(e  /El^]  (35) 


and 


UUn  ± (-  - l-l"')'  - (H  - 

E^  = — E (EE  )?.  (36) 

+ (M^  - 1)^ 


We  note  that  the  minimum  value  of  jj  after  two  scatterings  is  given  by 


(i  = U - 

m o 


(37) 


The  expression  for  can  now  be  written 


Y2(E,x,h)  = ~ U(x)U(u  - ^^) 


E 

X J‘  +dE'?(Ep,E')5(E',E)T2(^g,-t(E'),L(E,n);x) 


X [g(E  ,e')  A(eSe^  . 


(38) 
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D.  The  Heru'lLe  Polynomial  Expansion  of  ^(E,x,^) 

We  solve  Eq.  ( .')  by  expanding  ‘5(E,x,^)  in  llcrmiLe  polynomials  in 
the  variable  x and  Legendre  polynomials  in  llie  variable  Starting 

with  the  Hermite  expansion  we  write 


and 


$(E,x,  a)  = 


exp 


} . 
1 


X - x(E;\-  ( 

5(E)  ) ( 


X 


1 


X - x( E ) N 

5(E)  J 


(33) 


’i5.(E,x,u) 


e cp 


X 


('  ) 


where  h 

s 


is  the  normalizing  factor 


h = 
s 


/7i  2'"  sJ 


(M) 


arising  from  the  integral 


r d y e ^ 11  (y)  H (y)  = It  6 
t s n ssn 


(’t2) 


We  substitute  the  expansions  (39)  and  (UC)  into  Eq.  ('3),  use  the 
fact  that 


(liS) 


il 


\ 


and  operate  on  both  sides  with 


-C^ 


d X 


exp 


i:(E)  / 


and  obtain 

- ■ n - ' -(E) 


= ^ dE  ' J(E',E)  r " djj'  dyV'-' 

'"E  -1 

(...) 

X 5(  V • V ' - g(E  ^.E) ) 


LQ  (E',E)  TY,  (E',u')  + r (E' 
n , s - ; ; s s 


y )] 


s = c 


where 


g (E'.E)  .1,; 


d X 


I.  r- 


n\^  §(E) 


P) 


C^5) 


The  evaluation  of  ^(E^,E)  is  left  to  Appendix  B. 
Henceforth,  we  will  use  the  notation 

7 = x(E),  § = §(E),  x'  = x(e'),  5'  = P(E')  . 


(i+6) 


The  functions  x(E)  and  |(E)  are  yet  to  be  specified.  The  advantage  of 
introducing  x(E)  and  §(E)  is  that  much  of  the  gross  behavior  (as  a 
function  of  x)  of  the  distribution  function  can  be  built  into  the 
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Gaussian  factor  in  expansion  (39).  Consequently,  many  fewer  terms 
should  be  required  in  the  llermite  expansion  than  would  be  the  case  if 


x(E)  = 0,  §(K)  = constant  . 


E.  Tlie  Legendre  Polynomial  Expansion  of  5(E,x,  i,;) 

We  now  expand  ^ ~ S^E^E))  in 

Legendre  polynomials;  specifically, 


rg(E,a)  = U 


E2m  ■* 
2 


9 (E)P  (u)  , 
^s  ,m  m 


m = 0 


/r.  V V ' 2m  ^ 1 

V (E,u)  = u > 5 T ^ (E)P  (u) 

J ^ y,s,m  m 


m = 0 


6(0  . 0'  - g)  = ^ P (g)  [P  (^)P  (n') 

<=•  m m m 


m = C 


111 

^ £ (I'l  J)j 


r = 0 


By  substi  tutriig  these  expansions  into  Eq.  (4^1)  and  performing  angular 
integrations,  we  obtain  a coupled  set  of  integral  equations  for  Llie 


■''n.d"’  * “n,d")  ■'  .^''''<'E'S(E',E)P^(g<E',E))Q_,___(E%E).p^^,<E')  , 
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t 


whore 


V^^f(E)  = ^ ciE'^(E',E)P^(g(E',E)) 

s = 0 ^ 


(52) 


1 1 -1 

J,,^^(E)  = ^ J^^dE'‘/(E',E)P^(g(E% 


E)) 


s = 0 


X Q„  (E',E)y„  ,(E')  . 


n,  s 


(53) 


An  explicit  expression  for  Y _ ,(E)  is  derived  in  Appendix  C.  Although 

Eq.  (51)  describes  a set  of  coupled  equations,  the  coupling  is  such  that 

the  c,  .(E)  may  be  obtained,  order  by  order,  in  a natural  sequence.  The 
n ^ “C 

equations  for  the  functions  ^ „(E)  are  uncoupled^ 


?(E)Nj'(E)9  ,(E)  = V (E)  -t  JP^dE'y:(E',E)Q  (E',E)^  (E'),  (>4) 

o , c.  u ^ L/  w u ^ c ^ ; 1- 


and  so  these  functions  may  be  obtained  directly,  say  for  a sot  t = 0,  1, 

. . . , 1.  The  equations  for  tlie'+i  ,(E)  depend  on  9 ,m  = 0,  1,  ..., 

t 1,  and  so  y . may  now  be  found  for  the  sot  of  i values  T = 0,  1, 

1,-L 

.t  - 1.  This  process  is  tlien  continued  until  q/_  is  determined. 

f.,o 

If  at  that  point  we  see  that  more  terms  are  needed,  we  can  solve  suc- 
cessively for  > _ , y _,•••,  ■r_  without  modifying  the  solu- 

c,-f.  + 1 1,-f-  1+  1,0 

tions  for  the  terms  obtained  previously. 
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Kquation  ( is  modified,  of  course,  wlien  the  inelastic  loss  term 


in  the  Boltzmann  equation  is  taken  into  account.  The  form  of  the 
equation  remains  unchanged  however  and  the  procedure  for  its  solution 
is  unaffected.  Details  of  tlic  changes  introduced  by  including  the  con- 
tributions of  the  inelastic  losses  are  given  in  Appendix  A. 

Tlie  matrix  inversion  method  for  solving  the  integral  equations  for 

the  ■ , is  described  in  detail  in  Appendix  D. 

n,T 


F . Calculation  for  200  keV  ;\ntimony  Atoms  Incident  on  Germanium 

In  our  calculation  we  used  the  Winterbon  analytical  fit  to  tlie 

LNS^°  elastic  cross  section,  with  however  two  modifications.  The 

Winterbon  fit  to  the  elastic  cross  section  diverges  at  small  energy 

transfers  as  we  have  therefore  introduced  a low-energy  transfer 

cut  off  T,  at  lli.5  eV,  a reasonable  estimate  for  the  energy  needed  to 

displace  a germanium  atom.  With  this  cut  off  however,  the  total  cross 

section  goes  to  zero  precipitously  at  T . As  is  pointed  out  in 

Appendix  D,  this  introduces  difficulties  with  the  numerical  solution 

for  ..  We  therefore  introduced  one  furtlier  modification  in  the 

r\,l 

cross  section  below  100  eV  which  allows  the  total  cross  section  to  go 
to  zero  smoothly  with  energy,  iloth  modifications  are  described  in 
detail  in  Appendix  E.  We  use  the  LSS^’  evaluation  of  the  inelastic 
stopping  power,  but  again  with  a modification  introduced  for  numerical 
convenience  at  energies  below  100  eV.  These  alterations,  because  they 
are  made  at  energies  so  low  compared  to  the  incident  energy,  do  not 
alter  quantities  of  physical  interest  such  as  range  and  deposited 
energy  distribution  to  any  significant  extent. 
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As  stated  earlier  we  have  found  it  expedient  to  neglect  the  effect 
of  inelastic  energy  loss  on  anti  Y.  • The  depth  distribution  of 

electronic  energy  loss  near  tlie  surface  will  of  course  be  adversely 
affected,  but  Y itself,  and  therefore  elastic  losses,  is  hardly  affected 
at  all.  As  an  illustration,  in  Fig,  1 we  plot  Y^  integrated  over  energy 
(the  energy  dependence  is  set  by  the  6-function)  as  a function  of  angle 
lor  t!ic  case  where  inelastic  losses  are  included  and  also  \diere  they  are 
neglected.  Albeit  it  is  a log  plot,  the  results  arc  indistinguishable. 

The  functions  x(E)  and  F(E)  were  determined  by  a quasi-self- 
consistent  method.  A pliysically  reasonable  form  was  initially  chosen, 
and  Eq.  (jl)  (modified  according  to  Eq.  (All’))  was  solved  for  the 

.(El.  Moments  of  the  range  distribution  of  beam  particles  having 

n,l 

energy  E were  then  calculated  (see  Appendix  F)  and  a new  analytical 
fit  for  x(E)  and  F(E)  obtained.  No  furtiier  iteration  was  found 
necessary.  Figure  2 depicts  the  energy  dependence  of  x and  § found  in 
this  way. 

Expressions  for  the  distribution  of  stopped  ions  and  for  the  dis- 
tribution of  energy  deposited  in  elastic  collisions  are  derived  in 
Appendix  F.  It  is  shown  there  that  these  distributions  depend  only  on 
and  The  range  and  deposited  energy  distributions  sho\,m  in 

Figs.  3 and  h include  contributions  for  values  of  n up  tlirough  7-  At 
tliat  point  the  Hermite  expansion  has  not  yet  quite  converged,  but  the 
resultant  distributions  should  be  good  to  about  K'^. 
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LOG 


tNtRGY  iNTeGRATED  SINGLE  SCATTERED 
HUX  WITH  AND  WITHOUT 


INCLUSION  OF  INELASTIC  LOSSES 


O k 1364  LINDHARD  UNITS 


Kiu.  1 “ flux  of  once  scattered  particles  inte^ated  over  onorjjy  with 
mid  without  the  inclusion  of  inelastic  losses.  An  inelastic  loss  proportional 
to  velocity  is  used,  with  a value  of  k of  .1354  Lindhard  units. 
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Appendix  A.  Inclusion  of  Inelastic  Losses 


As  indicated  in  the  text,  inelastic  losses  (in  the  continuous 
slowing  down  approximation)  are  taken  into  account  hy  including  a term 

I(E,x,u)  = - (S(E)?-(E,x,u))  (Al) 

on  the  left  hand  side  of  Eq.  (23),  where  S(E)  is  the  inelastic  stopping 
power.  In  order  to  obtain  the  alterations  made  by  this  term  to  Eq,  ( .1) , 
we  expand  i in  Hermite  and  Legendre  polynomials.  First  we  make  the 
Legendre  expansion  sho^>m  in  Eq.  (^t8)  , 

UE,x,u)  = ^ q-^(E,x)P^(Li)  . (Ai:) 


We  divide  by  2/(21  t 1)  to  find  the  contribution  of  the  -C  th  harmonic  to 
the  transport  equation  is 

I^(E,x)  = - ^ (S(E)q-^(E,x))  . (A3) 

We  now  introduce  the  Hermite  expansion  of  , 

'i',(E,x)  = h ‘-y  ,(E)e'y  H (y)  , (A4) 

I ’ / j m m,-t  m 

m = C 


where 

y = (X  - x)/r  . 


(A5) 


Now  the  Hermite  polynomials  satisfy  the  relation 


d 

dy 


(H  (y)e  ^ ) 
m 


H 


m ♦ 1 


(y)e 


-y‘ 


(a6) 


and  therefore 
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I (E,x)  = - S'(E)  h (E)  H (y) 

/ j m •{  m 


m = C 


- s(E)  ii  , 

m = C 


' (E)  e'^  H (y)  • ^ 


' E 


h-V  ,(E)  e " II 


m ' " ' / j m , 1 m 

m = 0 


X 

F 


- / 


(X  - >^) 


Using  the  recurrence  relation 


= *‘m  + i^y)  " ‘^m-  .^y>  ^ 


we  obtain 


I^(E,.,  . =->’'  ^ {:  - S'(E)V^_^(F.)  - S(Ov;_^(El 


m = 0 


V /EiH  (y) 
c rn^-t/  m L 


— y (E)  H (y)  ] 


Operating  on  I.  by  ^(E)  ; dy  H (y),  we  find 
-t  «-  n 


, i^y)  ^ 


(A7) 


(A8) 


(A9'l 
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(AIG) 


e' 


^ f ^ (n  ' 1)-;^  /E)  + S(E)q>'  (E) 


x' 


Now 


h = /rf  2^n  I , 
n ’ 


thus 


* 2n(n  - l)-l^c|,^__^^^(E))  . 


(Alx) 


(A12) 
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Appendix  B.  Evaluation  of  ^(E^,E) 


In  this  appendix  we  evaluate  the  function 


'v.-y 


/2 


(Bl) 


which  appears  in  Eqs.  (^^4)  and  (51)  and  where 


h = /n  2®s  : 

S 


y = (x  - x)/5 
y'  = (x  - x')/§', 


(B2) 


and  where  the  are  Hermite  polynomials  of  order  n. 
VJe  first  rewrite  Eq.  (Bl) 


(B3) 


We  remark  that 


g(E  ^E)  -OyS>n, 


(Bl) 


and  use  the  fact  that 


2 = H__(x  , y)  . <»5> 


m = 0 


to  obtain 


2l 


I 

I 


I 


X r”dy'e^  H(y')H  (/cy' 
^ s m 

-00 


UiC) 


T!ie  integral  can  be  done  explicitly,''-  and  with  a cliange  oC  summation 
index  m = s 4 2r,  we  obtain 


Qn^3(E',E)  = (/2 


F2(— - 1]'^  X 


( n - s - 2 r ) .' 


r =0 


X H o - I - ) , 

n - s - 2r  r ’ 


(BO 


where  the  upper  limit  of  the  summation  [(n-s)/2"  denotes  the  greatest 
integer  less  than  or  equal  to  (n-s)/2. 
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Appendix  C.  Evaluation  of  : (E) 

liy  definition  i ,^^^(E)  is  related  to  Y (E,x,u)  by  the  expression 


= J d^  1 (U)§  J^dx  r^(E,x,a)  . (Cl) 


-1  r.  ® 


X - X . 


To  obtain  a numerically  useful  expression  for  1 . however,  the  best 

^:n,l  ’ 

method  is  not  the  most  direct  one  where  T,  is  taken  from  Eq.  (38). 
Rather  it  is  better  to  start  with  the  defining  equation  for  V^(E,x,^) 


^ T*  ^ ^ 

u •..  + ? (E)v^  = dE'y(E',E)v,(E',x)  X 

E 


X J dv'6(v  • v'  - g,)  6(  iu'-  g^,)  , 


(C2) 


where 


g.  = g(E',E)  and  g^  = g(E|^,E') 


(C3) 


and  g(E^,E^)  is  the  cosine  of  the  laboratory  scattering  angle  for  a beam 
atom  entering  with  energy  E^  and  leaving  with  energy  E^.  After  we  expand 
the  6 -function 


6(v-v'.gj= 


1=0 


' (t^  m(9-cf')' 


(C4) 


m = 1 


the  integrations  over  ' and  pi'  can  be  done  imr.iediately  to  yield  for  the 
RHS  of  Eq.  (C£) 
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( 

I 


‘i(K  ,K)Vi(E  ,x) 
RHS  = d E 


g(E, 


1<.E  ,x) 

,E')  ^ 


2-L+  i 


2 I 


P/gi)P/g.  )m  Pc  (u).  (C^) 


Now  expand  Y, , 


Y^(E,x,^)  = H ^ 


(c6) 


I f we  note  that 


a X 


1 5 „ /p 

2 3x 


i = C 


X)  u'^P^(u) 


(C7) 


^2  ax  '*'2;^  - 1 2 Sx  '2;t  t 1^  ’ 


^ = C 


we  obtain 


2t  • i ax  '^2;-f,  - 1 *■  2-t  < i ax  '^2;t  t 1 


3E,  , 3=(E',E)V  (E',x) 

r ‘IE'  


(C8) 


We  now  expand  Y (E,x)  in  Hermite  polynomials.  We  let 


y = (x  - x)/5  , 


(C9) 


and  write 


n = 0 ’ 


(CIO) 
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where 


h = /tt  & n'  . 
n 


We  substitute  this  sum  for  Y in  Eq.  (08),  note  tliat 

2 ' 


^ V=-»  “ - !<>'> 


and  use  the  orthononnal ity  property  of  the  Herraite  polynomials, 


;*dy  e-y  ii3(y)H^(y)  = 


to  obtain 


(-  2 s)  [ 


. 4^y 


2-t  1 i r-;s  - 1,  f,  - 1 2-L  I 1 2;s  - 1,  1- 1 


3E  5(E',E)P^,(g  )P  (g  ) a, 

= clE'  i-.i  ...-^  ^ r dxf  (E',x)H  (y) 

"E  g2  ® 


Taking  V from  Eqs.  (26)  and  (27)  and  noting  tliat 


- f>= 


k = 0 


and 


where  the  symbol  [k/2]  denotes  the  largest  integer  contained  in  k/2, 
find  that 


(Oil) 


(012) 


(013) 


(Oil) 


(C15) 


(016) 

we 
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d x;^(E',x)H^(y')  = 2 


s : 


7H_.(-/T|)  X 


(s  - k) ; s - k 


Z(-  1)^  _21 

V.' 


V = 0 


. ov,  k-2v  » u 
' k -^v  ^ ^ 

VTTF); 

u = 0 


(Cl7) 


where 


= a'^(Eg)  and  I(e')  = a^(E')/g, 


(C13) 


The  final  result  then  is 


(-2  s) 


I 


t + 1 


2(,+  l '^2;s  - l,f.  j ■*  ^ ^ ^ '^£;  s 


s; 


- t E E 

B k = 0 - V = 0 


(-  j.)  / 2^  . 

v;  Ki^zJ 


k - 2 V 


X (Cl9) 


k - 2 V 

e;' 

a = 0 


„ ?(E  ,E')5:(E',E) 

J r.  / <2 


d E 


a (E  ) 


We  have  therefore 


(E) 


1 odE  , 3'(Kr2E')?(E',E) 

tTET  J_  TTT;: P^(g(Ep,E'))P^(g(E',E)) 


(C20) 


5a^(E)tg  '^F.  O^(E') 


and 


la  (E)Y  (E) 
s\n,l 


(C21) 


= 2n 


[2^  ’sin  - 1,(  - 1<«  * sHl  - 1,(,  . l‘«  j ' \,l 
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Appendix  D.  Nimiorical  Solution  of  ihc  Integra]  Equation 

In  this  appendix  we  describe  how  the  integral  equation  for  y (E^ 

' }' 

is  solved.  This  will  serve  as  an  illustration  of  how  tlie  equations  for 
all  tlie  Y „(E)  are  solved.  For  the  purposes  of  this  section  we  will 

n j'u 

neglect  inelastic  losses  (but  will  indicate  at  the  end  how  they  are 

treatedl  and  will  assume  the  cross  section  is  given  by  the  Winterlion, 

Sigmund  and  Sanders  analytical  fit  to  the  Lindhard,  Nielsen  and 

Scliarff'  cross  section.  Since  the  LNS  cross  section  becomes  singular  as 

the  energy  transfer  goes  to  zero,  we  introduce  a cutoff  T to  represent 

a minimum  energy  transfer  allowed  in  an  elastic  collision.  Tlie  metliod 

described  below  works  equally  well,  moreover,  witli  otlier  cross  sections. 

With  the  use  of  the  LNS  cross  section,  the  function  v (E)  satis- 

0,0 

fies  the  equation 


^,(E1  = V (El 


t ^ d e'  U(E„  - 3 T 
-E  * T ® 


e')5'?(e',e')9^^^(e'i  , 


(Dl) 


where  all  of  the  terms  appearing  here  have  been  defined  in  the  body  of 
this  report.  T f we  define 

m(E)  = min  rp  E,  Ep  - 3 T]  , ( D£) 

then  we  have 

(D3) 

4 dE'E'‘J(E',ElV^  „(e')  . 

V LJ  j U 

E 1 T 
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Briefly,  given  an  appropriate  energy  mesli  {B.},  oiir  aim  is  to 

express  Eq.  (Dii')  as  a set  of  linear  equations  in  the  set  {y  (E.)], 

whicti  we  will  tlien  solve  by  straightforward  matrix  inversion. 

Consider  the  integral  in  Eq . (Do).  We  carry  out  the  integration 

using  Gauss-hegendre  quadrature.  Thus  if  {i:.}  and  fw.}  are  the  sets 

1 ^ 

of  Gauss-Legendre  ordinates  and  weights,  respectively,  for  n-point 
integration,  then  the  integral  in  Eq.  (D3)  may  be  expressed 
. m(  E 

E.  ‘ 

J 


(D'O 


where 


and 


= B(E.)  + R(E.)z. 

Ji  J J 1 


R(E)  = I [m(E)  - (E  I T)  J,  B(E)  = ^(m(E)  -t  E i T)  . 


(d:^) 


(Di'O 


This  method  fails  when  singularities  (or  near  singularities!  occur 
in  the  integrand;  when  such  is  the  case,  the  singularities  should  be 
removed  before  integrating.  For  example,  tlie  WSS^  representation  of 
the  LNS^°  cross  section  is  proportional  to  (e'  - tliat  is, 

■?(E',E)  = G(E',E)/(e'  - E)-^/-’  , (D,) 


where  G(E',E)  is  well  beliaved.  In  this  case  the  singularity  is  removed 
by  introducing  the  change  of  variable  y ^ (E^-E)”^^'\  so  that 
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-in  (El 

d e'J(EV?(F/,E)<(.  , (E') 

C I o 

E'T 

(M) 

= •5  d y 5(E  I y“'’)'ji(E  I y*'',E)t-  .(Ely"'). 

[iiKE^-  El- 


The  new  intogrnl  contains  no  singularities  and  can  be  evaluated  as 
described  above.  Another,  not  so  obvious,  singularity  is  associated 
with  the  LNS  cross  section  and  with  the  cutoff  T ; as  E approaches 
T/7,  where 


(1  I M)^ 


(D3) 


tlie  total  cross  section  c ( E)  rapidly  goes  to  zero.  The  loft  side  of 

Eq . (DMl  remains  finite,  however.  This  means  that  r (E)  is  increasing 
T 

just  as  fast  as  o (E)  is  decreasing.  We  have  handled  this  situation  by 

modifying  the  LNS  cross  section  at  energies  less  than  ICO  e\'  so  that  tliis 
T 

sharp  cutoff  of  cr  (E)  does  not  occur.  The  modification  is  described  in 
Appendix  E. 

Let  us  assume  that  all  singularities  liave  already  been  removed  in 
Eq . (DL).  Then  for  eac'.i  mesli  point  E.  we  have 


r(E  )n  (K  )r,  , (E  ) = V (E.) 
J J ' J J 


(fjP . 


(DIO) 


Generally,  the  argument  of  ^ on  the  right  side  will  lie  among  the 
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mesh  points  [•‘■j]-  cubic  spline  interpolation  to  express 


V,  "’ji'  ■ ,„<V' 


(Dll) 


where  the  coefficients  are  determined  by  the  cubic  spline  equations 
and  found  without  having  to  know  tlie  values  t (E,  ).  Equation  ( DIG)  , 

wliich  now  describes  a set  of  linear  equations  in  [G  (E,  )1,  therefore 

' ,c.  k ^ 

has  the  form 


(D12) 


where 


A..  = 5(E.)o  (E.)6.. 

jk  = j J jk 


Ttie  desired  solution  is  then 


(D13) 


9 „(E,  ) = y AT:  V (E.)  . 

c ,o  k ^ Tcj  c,o  y 


(Dll.) 


As  a test  of  the  accuracy  of  the  method,  we  considered  the  equation 


, m(E) 

D(E)  = V(E)  4 ^ d e' K(E',E)D(E')  , 
E+T 


(D15) 


wi  th 


K(E',E)  = [ 1 + (E'-  E) ]/(E'- 


(D16) 
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The  function  V(E)  is  a ratlier  complex  function  of  E and  T,  siicli  that  the 


solution  of  Eq . (Dli")  is 


D(E)  = (E|j  - o T - E)-  . 


(D17) 


By  employing  tlie  metliod  described  above,  a mesh  of  h2  energy  points,  and 
ZP-  Gauss-Legendre  points,  we  reproduced  the  exact  solution  at  all  the 
mesh  points  to  an  accuracy  of  better  than  per  cent. 

In  Appendix  A,  we  have  shown  that  inelastic  losses  are  taken  into 
account  by  adding  a term  I ,(E.)  to  the  left  hand  side  of  Eq.  ( Dl(  ) > 
where 


I ,(E)  = - E(E)[[S'(E)  + (n  t 1) 


S (E)]cf  (E)  + S(E)cp'  (E) 
n,f.  n,/. 


— / 

+ 2n  ^ S(E)cf  (E)  + 2n(n-l) 


y 


(D13) 


where  S(E)  is  the  inelastic  stopping  power  and  the  "prime"  denotes  the 

derivative  with  respect  to  energy.  From  the  cubic  spline  equations,  we 

obtain  ' . ( E .) 

n,-f.  J 


(E.) 


L 


Ji  n,l  1 


(D19) 


The  resulting  equations  can  then  be  solved  as  indicated  above  by  straight- 
forward matrix  inversion. 
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Appi'inlix  K.  Modification  of  LNS  Cross  Section  and  LSS  Inelastic 
Stopping  Power  at  Low  KnerKV 

. Cross  Section  Modification 

The  uniiioditied  LNS  elastic  scattering  cross  section  is  given  by 


f (t^/^)U(PE  - E')  , 
LNS  27  E,'"  ' 


(El) 


where  a is  tlie  Lindhard  screening  length,  E^  is  the  Lindhard  energy, 


^ ,T  = e'-E, 


(E2) 


f(t'^  1 is  given,  in  the  Winterbon  approximation,  by 


f(t^/')  = 1.7^9  t’/'-'  [ 1 + (2.610  t'  , (E3) 

tlie  step  function  U(PE-E^)  merely  expresses  the  kinematical  elastic 
scattering  condition  that  the  initial  energy  of  the  scattered  particle, 

E,  is  greater  than  the  minimum  possible  value  of  the  outgoing  energy. 

All  other  symbols  are  defined  i i the  body  of  the  paper.  The  cutoff  at 
low  energy  transfers  is  introduced  by  multiplying  Eq.  (El)  by  another 
step  function 

>’  (E',E)  = - E'f  'V  f(tV“'’)U(pE-  E'1U(E'-  E - T ) , (El) 

LNS  2 7 Ej 

where  is  the  cutoff  energy  transfer.  For  expediency  in  the  numerical 
solution  of  the  integral  equations,  we  also  wanted  the  total  cross 
section  to  vanish  smoothly  at  low  energies  and  so  wo  introduced  one 
other  modification  for  energies  less  than  an  energy  E (which  we  took 


3b 


J 


to  be  IOC  eV) 


. / EE 

U(E'  - r ) 

? Cf'  F'l  E - T 

U(E'-  E - T ) 

c 


(E^) 


, E' < E , 


In  Fig.  J we  show  the  energy  dependence  at  low  energies  of  both  the 
modifleo  and  unmodified  LNS  total  cross  sections  (both  including  the 
low  energy  transfer  cutoff  at  eV) . 

Inelastic  Stopping  Power 

The  LSS  elastic  stopping  power  is  proportional  to  the  particle 
veloci ty 

S(E)  = KE'-''^’, 


(e6) 


where,  for  the  case  of  antimony  on  germanium,  we  take  K =1.32  MeV^^^/^m. 
For  numerical  convenience  we  have  added  an  additional  term  at  low  energy 


(E() 


so  that 


S(E)  = 


K E- / , E > E 

I E , 


(E8) 


where  again  we  have  taken  E -=  1<X  eV. 
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Appendix  F.  Range  and  Deposited  Energy  nistribut  ions  and  Tlieir  MomeiUs 

1.  Distribution  of  Stopped  Beam  Ions 

As  stated  in  the  body  of  the  text,  we  have  written  Y as  tiic  sum  of 
four  terms. 


Y = Y^  * Y,  ^ Y.  It,  (FI) 

representing,  respectively,  the  flux  of  unscattered,  once  scattered, 
twice  scattered  ions  and  the  flux  of  ions  scattered  three  or  more  times. 
Because  the  energy  transferred  in  an  elastic  collision  by  an  ion  of 
energy  E can  be  no  greater  than  /E  , no  contribution  to  the  distribution 
of  stopped  ions  can  come  from  Y^.,  Y^^,  or  Y.^.  Con;ider,  therefore, 

Eq,  (23),  the  equation  for  9 (and  include  the  inelastic  contribution) 

u + Na^(E)i|  - (S(E)§) 


(F.) 


,PE 


2tt 


^E 


dE'5f(E',E)  r dC'6(v-v'-  g(E',E))u'’^  [ Y.(E',x,^^  $(E',x,u'", 


The  X component  of  the  current  of  particles  scattered  three  or  more 
times 

Er-ST 

dE$(E,x,n) 


j(x)  = 


-1 


must  obey  a particle  conservation  rule 


^ J (x)  do  (x)  d 1 (x) 

— + ^ 

d X d X d X 


(E'O 


where  do(x)/dx  is  the  distribution  of  stopped  beam  ions  and  dr(x)/dx 
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is  the  distribution  of  twice  scattered  particles  wtiich  are  elastically 
scattered  a third  time  in  the  interval  between  x and  x ♦ dx  . I f we 
operate  on  Kq.  (F.t)  wi  tii 


\ 

J 


1 

d Li 

1 


-1 


Eb-3T 
> d E 


make  use  of  tlie  relation 


kTT 


TT  . . 

d¥'6(vv'-  g(E',E)) 


= 1 , 


(F^,) 


and  recall  'he  definition  of  the  total  cross  section  from  Eq.  (15) 


No^(E) 


, E 

\ d e'5:(E,E')  , 
E/p 


we  find 


(f6) 


^ j (x) 

d X 


Eb-3T 
. - 1 ' . _ 


^ dL.U  ^ d E i^(S(El^(E,x,^)) 


(FT) 

^ Eb-3T  PE  1 

dE  \ d e'  \ duV'-' '?(E^E)Y^(E^x,^')  = 0 . 

w s) 

^ E -1 


We  find  therefore  that 


d<"  ( x) 

d X 


Eg-iT 


(F8) 


We  expand  $ in  Legendre  polynomials 
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^(E,x,u') 


(F9) 


and  find  that 


Eb-3T 

d E -2-  (S(E)c?  (E,x)) 
SE  0 


= S(o)9  (o,x) 
o 


(FIO) 


since  y vanishes  for  E = E - 3 T.  The  Hermite  polynomial  expansion 

O 

gives 


dg  ( x) 

d X 


r 

S (c  ) exp  1^- 


(Fll) 


1 X - x(r  ) 

V §(o)  ) 


] 


h--  H 
n n 


X.:  x(ol 

§(c)  ) 


V ,(c) 
n,o 


It  is  clear  that  without  inelastic  losses,  there  would  be  no  stopped 
ions  - any  ion  whose  energy  falls  below  T would  have  no  mechanism  for 
further  energy  loss  and  would  never  come  to  a stop. 

Moments  of  the  range  distributions  are  easily  calculated, 


f>  , m do  (x) 
R = d X x — r-^ 
m J dx 


= S(ol  y C (o)y  (ol 


(F12) 


where 


C 

n,n. 


(E) 


h*  ^ d X X 
n 


H 

n 


..-.xCEJ 

5(E1 


) exp 


X - xCEl  \ 
§(E)  1 


(F13) 
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UsinK  the  recurrence  relations  for  llermite  polynomials 
2yll^(y)  = ' 2nH^_  ^(y)  , 


and  recalling  that 


h = /tt  £ nJ  , 
n ’ 


we  arrive  at  a recurrence  relation  for  the  C . I f we  let 

n,m 


C (E)  = §(E)  C (El  , 
n ^ m n ^ TT1 


then 


C = x(E)  C , I (n+l)5:(E)  C , (E)  +|?(E)(l-6  )C  , 

n,m  n,m-l  - n + l,m  + l'  c;  n-1  , 


I f we  note  that 


C =1  and  C = 0,  m<  n 
o > o n , m 


we  can  calculate  all  the  C (E). 

n,m 


2.  Energy  Deposition 

The  energy  current  T|(x)  is  given  by 


11(x)  = d E EY(E,x,n) 


-1  o 


if  we  define  X(E,x)  hy 


X(E,x)  - ^ d^L  Y(E,x,^)  , 


(Fl^) 


(Fl^) 


(Fl6) 


(E) 

'1  • 

(FIT) 


(F18) 


(F19) 


(F20) 
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t lu'M  till'  B<^itzmann  equation  yields,  after  a little  manipulation 


dTl(x)  ’■-B 

— + d K T(E)X(E,x)  I \ d E S(E)x(E,x)  - E„6(x)  , 

•J  ^ O 

C 


T(E)  = d e'(E-  e')7(E,e')  . 
E/3 


We  tluis  can  identify  the  quantity 


d X 


\dECT(E)  I S(E))X(E,x) 
o 


(Fi-1) 


(F22) 


(F2o) 


as  the  depth  rate  of  energy  deposition.  The  contribution  from  elastic 
and  inelastic  collisions  are  clearly  separated.  Using  Eq . (10)  from 
the  n'.'.in  body  of  the  text,  we  have 


deixl 


Eb 


o ^ 


^ r~  r*  — 2*'  00  

+ \ d E[  T(E')  I S(E)]exp[^-  ( " ) i £ h'^l  i (E) 

J L _i  L V i(t)  / J n n \ §(E)  /L  n,o^ 


+ V (E)  I 
P?;n,c^  J 


(F^l) 


The  moments  of  the  energy  distribution  are  given  by 

,,  s d£(x) 

V = \ d X X — 

s J dx 


E,/P 


£ Vd  E TlfE)  •»  S(E)lc  (EiT-.  (E)  + 9 (E)  I, 

n = 0 J J n,s  L j';n,o'  ^n,o'  'j’ 

o 


(F2^) 


(F2'6) 
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